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An  empirical  microscopic  recombination  model  is  developed  for  the  direct  simulation  Monte  Carlo 
method  that  complements  the  extended  weak  vibrational  bias  model  of  dissociation.  The  model 
maintains  the  correct  equilibrium  reaction  constant  in  a  wide  range  of  temperatures  by  using 
the  collision  theory  to  enforce  the  number  of  recombination  events.  It  also  strictly  follows  the 
detailed  balance  requirement  for  equilibrium  gas.  The  model  and  its  implementation  are  verified 
with  oxygen  and  nitrogen  heat  bath  relaxation  and  compared  with  available  experimental  data 
on  atomic  oxygen  recombination  in  argon  and  molecular  nitrogen.  Published  by  AIP  Publishing. 
[http://dx.doi.org/10.1063/L4986529] 


I.  INTRODUCTION 

The  direct  simulation  Monte  Carlo  (DSMC)  method,  pio¬ 
neered  by  Bird  over  fifty  years  ago,1  for  most  of  that  time 
has  been  considered  the  main  numerical  tool  for  modeling  rar¬ 
efied  chemically  reacting  gas  flows.2,1  A  number  of  chemical 
reaction  models  have  been  proposed  over  the  years,  from  sim¬ 
ple  empirical  ones  such  as  the  total  collision  energy  (TCE) 
model4  and  the  quantum  kinetic  (QK)  model5  to  more  com¬ 
plex  approaches  that  include  vibration  coupling  effects,6-1 1 
and  to  state-to-state  and  other  sophisticated  methods.1213,16 
Due  to  scarce  cross  section  or  any  other  microscopic  or  non- 
equilibrium  reaction  data  at  the  time,  the  main  reference  base 
for  model  validation  used  to  be  temperature  dependent  reac¬ 
tion  rates.  A  model  was  usually  developed  from  or  checked 
versus  available  rates  written  in  the  Arrhenius  format.14  The 
primary  focus  for  the  model  development  so  far  has  been  dis¬ 
sociation  and  exchange  reactions  for  diatomic,  and  sometimes 
for  polyatomic,  molecules.  The  recombination  reaction  was 
largely  avoided  in  the  DSMC  computations,  which  to  a  signif¬ 
icant  degree  is  related  to  the  area  of  application  of  the  DSMC 
method. 

This  area  is  traditionally  limited  by  the  high  computational 
cost  of  the  DSMC  method  when  applied  to  low  Knudsen  num¬ 
ber  flows,  especially  when  these  flows  are  three-dimensional. 1 5 
Even  for  two-dimensional  flows,  modeling  is  usually  con¬ 
ducted  for  Knudsen  numbers  Kn  >0.001,  where  the  impact 
of  recombination  reactions  is  almost  always  minor,  so  that 
they  may  be  either  completely  disregarded  or  modeled  with 
a  simple  and  not  necessarily  accurate  model.  In  the  last  few 
years,  however,  the  rarefied  gas  dynamics  community  has  seen 
the  development  of  efficient  algorithms  for  modern  computer 
architectures16-1'1  which  dramatically  expand  the  area  of  prac¬ 
tical  applicability  of  the  DSMC  method.  These  methods  have 
already  been  tested  for  reacting  air  flows.20  Today,  modeling 
of  gas  flows  at  Knudsen  numbers  as  low  as  10-4,  and  possibly 
even  lower,  appears  within  reach. 
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Coincidentally,  the  development  of  high  performance 
DSMC  algorithms  occurs  at  the  time  when  significant  push  has 
been  made  towards  increased  accuracy  of  modeling  of  colli¬ 
sion  processes  in  reacting  air  (see,  for  example.  Refs.  21-25). 
These  studies,  which  mostly  use  molecular  dynamics  and 
quasiclassical  trajectory  (QCT)  calculations,  have  provided 
detailed  translational  energy  and  internal  level  dependent  cross 
sections,  as  well  as  non-equilibrium  reaction  rates,  for  oxygen 
and  nitrogen  dissociation  and  exchange  reactions.  Detailed 
QCT  results  for  the  specific  reaction  and  known  potential 
energy  surface  of  interest  always  provide  the  highest  accu¬ 
racy  when  available  and  when  the  computational  penalty  of 
implementation  can  be  overcome.  As  such,  the  cross-sectional 
information  may  be  used  directly  in  DSMC  simulations,12,13,16 
although  full  scale  translation-rotation-vibration  energy  tables 
would  be  prohibitively  large  and  thus  require  some  reduced 
models  such  as  Ref.  26.  Moreover,  multi-dimensional  hyper¬ 
sonic  flow  simulations  still  need  fast  and  simplified  mod¬ 
els  that  capture  the  key  physics  and  that  can  be  used  on 
reactions  for  which  detailed  QCT  results  are  not  available. 
In  this  case,  newly  available  reaction  data  provide  exten¬ 
sive  basis  for  validation  and  analysis  of  empirical  DSMC 
models. 

An  example  of  such  a  validation  effort  is  the  recent 
work,27  where  the  extensive  and  high-quality  recent  qua¬ 
siclassical  trajectory  calculation  work  from  several  groups 
on  nitrogen  molecule  collisions  with  N28,2'1  and  ISE, 16,23 
covering  thermal  non-elastic  and  dissociative  processes,  is 
used  for  DSMC  model  verification  and  parameter  adjust¬ 
ment.  Most  widely  used  in  the  DSMC  community  models 
of  dissociation  reaction  were  considered.27  The  main  conclu¬ 
sion  was  that  the  extended  weak  vibrational  bias  model8,41 
(called  the  bias  model  hereafter)  provides  the  best  agree¬ 
ment  with  benchmark  vibrationally  specific  dissociation  rates, 
while  significant  differences  were  observed  for  the  other 
models.  The  bias  model  requires  no  a  priori  assumption 
of  an  Arrhenius  rate  coefficient,  while  providing  realistic 
dissociation  reaction  rates  in  a  wide  range  of  temperatures.8 
It  reproduces  vibration-dissociation  coupling,  equilibrium, 
and  nonequilibrium  behavior  as  a  natural  consequence  of 


1 070-6631/201 7/29(6)70671 06/1 1  /$30.00 


29,  067106-1 


Published  by  AIP  Publishing. 


067106-2  S.  Gimelshein  and  I.  Wysong 


Phys.  Fluids  29,  067106  (2017) 


its  physics-based  cross-sectional  function  and  just  a  simple 
chemistry-based  assessment  of  whether  the  reaction  will 
be  strongly  favored  or  not.  Our  most  recent  analysis  of 
nitrogen  dissociation27  has  shown  that  the  bias  model  pro¬ 
vides  good  to  excellent  agreement  with  QCT  non-equilibrium 
rates  for  oxygen  dissociation  reactions  as  well.  The  obvi¬ 
ous  advantage  of  the  bias  model  is  its  ability  to  capture  the 
vibration-dissociation  coupling  process,  being  at  the  same 
time  easy  to  implement  and  numerically  efficient.  More¬ 
over,  it  has  the  potential  to  be  applied  to  reactive  processes 
where  no  reliable  non-equilibrium  or  equilibrium  rates  are 
known. 

One  significant  disadvantage  of  the  bias  model  is  the  lack 
of  a  compatible  recombination  model  since  the  extension30 
to  the  weak  vibrational  bias  model8  has  been  provided  only 
for  the  dissociation,  but  not  for  recombination,  reactions.  A 
compatible  model  in  this  case  is  a  model  that  at  minimum 
satisfies  the  principle  of  detailed  balance  for  collisions  in 
equilibrium  gas,  and  at  maximum  also  captures  the  corre¬ 
sponding  equilibrium  constant.  This  lack  significantly  limits 
the  applicability  of  the  bias  model  to  near-continuum  flows 
where  the  recombination  reactions  may  be  important,  makes  it 
difficult  to  provide  comprehensive  model  validation  and  com¬ 
parison  with  continuum  approaches,  and  complicates  its  use 
in  hybrid  simulations.  The  recombination  model  that  appears 
closest  in  terms  of  its  compatibility  to  the  bias  model  is  the 
approach8  developed  for  the  original  weak  bias  dissociation 
model.  But  even  aside  from  the  fact  that  the  original8  and 
extended10  bias  dissociation  models  differ,  and  a  straightfor¬ 
ward  extension  of  the  approach8  to  the  extended  bias  model 
does  not  seem  possible,  the  recombination  strategy  of  Ref.  8 
has  several  inherent  deficiencies.  These  are  the  violation  of  the 
detailed  balance  requirement  for  the  rotational  mode,  deviation 
from  known  equilibrium  rate  constants,  and  extremely  difficult 
implementation. 

Another  possible  route  for  the  addition  of  the  recom¬ 
bination  to  the  bias  model  is  a  robust  approach  proposed 
in  Ref.  31.  The  approach 31  provides  the  way  to  sample 
post-recombination  vibrational  energies  of  newly  created 
molecules.  Its  main  idea  is  to  record  and  tabulate  vibrational 
energy  distributions  of  dissociating  molecules  in  equilibrium 
gas  for  a  number  of  reference  temperatures,  so  that  for  each  ref¬ 
erence  gas  temperature,  there  is  a  corresponding  distribution 
function  of  vibrational  levels  of  dissociating — and  therefore 
recombining — molecules.  The  tabulation  is  conducted  in  a 
thermal  bath  prior  to  the  DSMC  flow  modeling.  Then,  during 
the  actual  flow  modeling,  the  pre-computed  tables  are  used  on 
a  per-cell  basis:  each  cell  has  its  own  temperature  and  thus 
uses  its  own  part  of  the  pre-computed  table  of  vibrational  dis¬ 
tribution  for  recombining  molecules.  Although  the  approach 
is  developed  for  the  TCE  model4  of  chemical  reactions  and 
the  Larsen-Borgnakke  (LB)  modek12  of  internal  energy  trans¬ 
fer,  it  is  likely  extendable  to  other  collision  models.  Inherent 
limitation  of  the  approach  is  related  to  its  reliance  on  equi¬ 
librium  sampling.  It  may  require  prohibitively  large  tables 
in  non-equilibrium  conditions  of  different  translational  and 
internal  temperatures,  especially  for  quasisteady  state  condi¬ 
tions2 1  with  significantly  underpopulated  tail  of  the  vibrational 
distribution. 


The  few  other  models7,11'14  for  recombination  reactions 
in  particle  methods  cannot  be  extended  to  complement  the 
bias  dissociation  model,  primarily  due  to  the  coupling  of 
the  vibrational  mode  with  its  dissociation  probability.  These 
models  were  developed  for  some  specific  dissociation  models 
and  cannot  provide  the  detailed  balance  to  the  dissociation- 
recombination  process  when  the  bias  model  is  used.  For  exam¬ 
ple,  the  early  recombination  models33,34  were  constructed  for 
the  TCE  model4  and  used  a  simplified  approach  to  the  after¬ 
reaction  energy  redistribution,  with  no  account  for  vibrational 
favoring  in  dissociation.  The  vibrationally  favored  dissocia¬ 
tion  (VFD)  model  takes  such  a  favoring  into  account  but 
provides  only  a  schematic  description  of  recombination  mod¬ 
eling  and  does  not  force  the  detailed  balance.  In  addition,  it  is 
only  suited  to  the  corresponding  dissociation  model,'  which 
has  been  shown27,30  to  be  less  robust  and  reliable  than  the  bias 
model.  The  more  recent  recombination  modek13  was  designed 
in  the  quantum  kinetic  framework,  and  as  such  does  not  include 
vibration-dissociation  favoring. 

The  main  objective  of  this  work  is  the  development  of 
a  recombination  model  that  complements  the  bias  dissocia¬ 
tion  model  and  satisfies  two  key  conditions:  the  number  of 
recombination  reactions  at  any  given  temperature  is  based 
on  the  correct  recombination  reaction  rate,  and  the  princi¬ 
ple  of  detailed  balance  at  equilibrium  is  taken  into  account. 
The  first  condition  means  that  the  recombination-dissociation 
process  captures  the  known  temperature  dependence  of  the 
corresponding  equilibrium  constant,  while  the  second  con¬ 
dition  implies  that  at  equilibrium  the  number  of  molecules 
dissociating  from  any  rovibrational  state  (/,  v)  is  equal  to  the 
number  of  molecules  recombining  to  that  state. 

The  proposed  approach  to  modeling  the  recombination 
process  is  straightforward  to  implement  and  may  easily  be 
adapted  for  most  other  DSMC  dissociation  models,  such  as 
quantum  kinetic  or  total  collision  energy.  Note  that  the  present 
work,  and  all  the  related  models  referenced  that  consider 
simple  5-species  air  chemistry,  operates  within  the  simpli¬ 
fying  assumption  that  considers  only  the  ground  electronic 
states  of  the  molecular  species.  It  is  clear  that  in  a  more 
general  case,  particularly  with  ionization,  excited  electronic 
states  will  play  an  important  role.  Excited  electronic  states 
will  greatly  add  to  the  complexity  of  dominant  reaction  mech¬ 
anisms,  including  inverse  predissociation  and  vibration-to- 
electronic  energy  transfer  paths  (see,  for  example,  Ref.  36). 
The  models  mentioned  in  this  work  (whether  implemented 
in  DSMC  or  computational  fluid  dynamics,  CFD)  would 
need  significant  revisions  in  order  to  be  considered  for  a 
flow  where  the  latter  processes  are  significant.  However,  the 
present  development  may  be  regarded  as  a  step  toward  bridg¬ 
ing  the  gap  between  the  most  accurate  continuum  CFD  models 
that  conventionally  use  both  recombination  and  vibrationally 
favored  dissociation  processes  in  modeling  hypersonic  flows, 
and  their  kinetic  counterparts,  which  often  disregarded  the 
recombination. 

II.  THE  BIAS  DISSOCIATION  MODEL 

The  key  objective  of  any  DSMC  chemistry  model  is  to  pro¬ 
vide  physically  adequate  number  of  reactions  in  the  range  of 
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temperatures  of  interest  and,  ideally,  use  reaction  cross  sections 
that  provide  good  agreement  with  available  experimental  and 
theoretical  data.  Although  it  is  possible  to  construct  a  temper¬ 
ature  dependent  algorithm  for  reaction  modeling,37  it  appears 
natural  for  a  particle  approach,  such  as  the  DSMC  method,  to 
use  the  available  information  of  translational  and  internal  ener¬ 
gies  of  colliding  particles.  This  not  only  improves  the  physical 
realism  of  a  model  but  also  provides  the  opportunity  to  per¬ 
form  detailed  microscopic  verification  and  validation  of  the 
entire  reaction  process.  The  vast  majority  of  the  DSMC  reac¬ 
tion  models  are  therefore  cross  section  and  energy  based,"  and 
model  the  bimolecular  reaction  of  dissociation  of  a  molecule 
A\A2  consisting  of  atoms  A\  and  At  and  colliding  with  a 
particle  X, 

A\A2  +  X  — >  A]  +  A2  +  X ,  (1) 


as  a  three-step  process.  First,  upon  selecting  the  collision  pair 
A1A2  and  X,  it  is  checked  whether  the  available  collision 
energy  Ec  is  large  enough  to  exceed  the  reaction  threshold 
Ed-  Generally,  the  definition  of  Ec  depends  on  the  chosen  dis¬ 
sociation  model.  In  the  original  weak  vibrational  bias  model,8 
it  is  assumed  that  Ec  is  the  sum  of  the  relative  translational 
energy  of  the  colliding  pair,  Et,  and  the  vibrational  energy  of 
the  dissociating  molecule  A1A2,  Ev,  that  is  E,  =  E,  +  E„.  In 
the  extended  bias  model,30  Ec  =  Et  +  Ev  +  ,  where  E,  is  the 

rotational  energy  of  A1A2. 

If  the  available  energy  Ec  is  sufficient  for  dissociation, 
then  the  algorithm  proceeds  with  its  second  step,  which  is 
the  reaction  probability  Pd  calculation  and  the  acceptance  or 
rejection  of  the  colliding  pair.  The  colliding  pair  is  accepted  for 
dissociation  with  a  probability  Pd  which  is  always  a  function 
of  the  energies  of  the  colliding  pair.  The  actual  form  of  the 
reaction  probability  is  determined  by  the  model  used.  In  the 
weak  vibrational  bias  model,8  it  is 


Pd  =  Ad 


Ed -E„\ 
E,  ) 


while  in  the  extended  bias  model,30 


Pd  =  Ad 


L  _  Ed  -  Ev 
\  E,  +Er 


(2) 


(3) 


Here,  Ad  is  the  calibration  constant  chosen  of  order  unity  to 
match  the  known  reaction  rate  constant  of  the  given  reaction, 
and  A,  is  the  parameter  controlling  the  degree  of  vibrational 
favoring.  Note  that  the  reaction  probabilities  Pd  are  instan¬ 
taneous  and  depend  on  the  particular  relative  translational 
energy  E,  of  the  colliding  pair  A  \A2-X  and  the  rotational 
and  vibrational  energies  of  molecule  A1A2.  Note  also  that  the 
extended  bias  model  provides  better  fit  to  known  equilibrium 
and  nonequilibrium  reaction  rates  than  the  weak  vibrational 
bias  model8  and  other  widely  used  DSMC  models  such  as 
TCE  and  QK.27  It  was  therefore  chosen  in  this  work  as  the 
baseline  dissociation  model. 

The  third,  and  final,  step  of  the  dissociation  model 
is  the  transformation  of  the  molecule  A1A2  into  atoms 
A 1  and  A2  and  the  energy  redistribution  over  the  reaction 
products.  The  energy  redistribution  in  the  bias  model  con¬ 
sists  of  the  removal  of  Ed  from  the  available  energy  Ec, 
and  the  sampling  of  the  newly  created  atoms  and  third 


particle  velocities  using  the  Larsen-Borgnakke  mechanism  '2 
and  the  considerations  of  the  energy  and  momentum 
conservation. 


III.  GENERAL  PROCEDURE  OF  THREE-BODY 
RECOMBINATION 

The  inclusion  in  the  reaction  probability  of  the  bias  model 
of  the  rotational  energy  of  the  dissociating  molecule,  described 
in  Sec.  II,  is  necessary  from  the  model  accuracy  perspective. 
However,  it  makes  the  recombination  model8  not  suitable  since 
the  equations  described  in  Sec.  II  E  of  Ref.  8  can  no  longer 
be  used  or  extended  in  an  analytic  form.  For  example,  the 
dependence  of  the  dissociation  cross  section  of  the  model8  only 
on  the  vibrational  and  relative  translational  energies  allows 
one  to  express  the  cross  section  cry;  of  a  recombination  into 
a  vibrational  level  i  from  an  initial  free  state  /  in  a  simple 
form 

o-fi(e)  =  cr  ei((e))Pjf, 

where  ( e )  is  the  symmetrized  relative  translational  energy,8 
crei  is  the  elastic  collision  probability,  and  P,y  is  the  dissoci¬ 
ation  probability.  While  the  approach  has  to  rely  on  energy 
symmetrization,38  it  still  provides  an  explicit  expression  for 
the  stabilization  of  an  orbiting  pair  of  atoms  into  a  specific 
vibrational  state  i. 

Such  an  approach  is  not  suitable  in  case  the  dissociation 
cross  section  also  depends  on  the  rotational  state  of  dissoci¬ 
ating  molecules  since  both  rotational  and  vibrational  states  of 
a  molecule  created  after  recombination  need  to  be  specified. 
It  is  however  possible  to  use  an  efficient  and  straightfor¬ 
ward  procedure,  described  below,  that  allows  one  to  satisfy 
the  detailed  balance  requirement  and  at  the  same  time  cap¬ 
ture  the  temperature  dependence  of  the  equilibrium  constant, 
which  was  not  enforced  in  Ref.  8.  In  this  procedure,  the 
process  of  recombination  of  two  atoms,  ,4 1  and  A2,  into  a 
molecule  A1A2  is  modeled  as  a  two-step  process  of  colli- 
sional  stabilization  of  an  orbiting  complex  (A1A2)  by  a  third 
particle,  X, 


Ai  +A2  ->  (AiA2), 

(4) 

(A1A2)  +  X  — »  A1A2. 

(5) 

The  model  is  built  to  satisfy  two  important  constraints. 
First,  the  number  of  recombination  reactions  at  any  given 
temperature  needs  to  be  based  on  the  correct  recombination 
reaction  rate,  which  is  determined  as  the  ratio  of  the  corre¬ 
sponding  dissociation  rate  of  the  bias  model  to  the  equilibrium 
constant  specific  to  the  dissociating/recombining  molecular 
species.  Second,  the  model  must  take  into  account  the  prin¬ 
ciple  of  detailed  balance,  which  states  that  at  equilibrium  the 
forward  rate  for  each  elementary  process  is  equal  to  the  reverse 
rate  of  that  process.  Applied  to  the  bias  model,  that  principle 
implies  that  the  number  of  molecules  dissociating  from  any 
ro vibrational  state  (./,  u)  is  equal  to  the  number  of  molecules 
recombining  to  that  state. 

The  proposed  recombination  algorithm  may  be  used  as 
an  addition  to  any  conventional  DSMC  collision  scheme,  such 
as  the  majorant  frequency  scheme39  or  the  no  time  counter 
scheme.2  It  is  added  to  the  cell-based  atom-atom  elastic 
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collision  algorithm,  and  the  complete  procedure  that  includes 
the  conventional  elastic  collision  and  recombination  reaction 
related  steps  is  as  follows. 

1 .  The  first  step  is  the  conventional  DSMC  elastic  collision 
procedure  for  atom-atom  interactions.  In  this  procedure, 
two  atoms  are  selected  randomly  from  all  atoms  located 
in  that  cell.  The  number  of  pair  selections  Ns  is  deter¬ 
mined  by  the  collision  frequency  in  the  given  cell,  which 
in  turn  depends  on  the  atom-atom  intermolecular  poten¬ 
tial  as  well  as  the  time  step  A t  and  the  number  density  n 
and  the  number  of  simulated  atoms  N  of  colliding  atomic 
species  A  i  andA2, 

Ns  =  siiAiNa2  (crT  g)max  Af . 

Here,  e  =  0.5  for  Ai  =At  and  1  otherwise,  crT  is  the  total 
collision  cross  section,  and  g  is  the  relative  collision 
velocity.  The  subscript  max  denotes  the  maximum  value 
for  the  given  cell.  The  selected  pair  of  atoms  is  accepted 
for  collision  or  rejected  based  on  its  relative  collision 
velocity  gAaAl  and  the  selected  intermolecular  poten¬ 
tial,  such  as  variable  hard  sphere  (VHS)  or  variable  soft 
sphere  (VSS).2  S  Typically,  the  acceptance  probability  is 
_  (crTg)AlA2 

a  ~  (o-Tg)wax  ■ 

2.  If  atoms  A\  and  A2  are  accepted  for  physical  collisions, 
then  either  an  elastic  collision  or  a  recombination  reac¬ 
tion  is  possible.  The  probability  of  recombination  is  Pr, 
and  the  probability  of  elastic  collision  is  1  -  Pr.  Here,  P, 
is  a  function  of  the  relative  collision  energy  Et  r(A\ ,  A2)  of 
pair  A] ,  A2,  computed  to  match  the  equilibrium  constant 
as  discussed  below.  Note  that  only  two  outcomes  are  con¬ 
sidered  in  this  work  for  atom-atom  interactions,  an  elastic 
collision  and  a  recombination.  A  selection  method  such 
that  of  Ref.  40  may  be  used  if  other  reactive  mechanisms 
need  to  be  included. 

3.  If  £%  <  P, ,  where  S/?,  is  the  random  number  uniformly 
distributed  in  the  interval  (0,  1),  then  there  should  be 
a  recombination  reaction  between  atoms  A 1  and  A2 .  In 
this  case,  the  internal  (rotational  J  and  vibrational  u ) 
energy  states  of  the  newly  created  molecule  A1A2  are 
chosen  from  the  list  of  states  of  recently  dissociated 
molecules.  The  internal  states  of  recombining  molecules 
are  assigned  sequentially  from  the  local  list,  as  described 
in  detail  in  Sec.  VI.  Importantly,  such  a  procedure,  along 
with  the  balanced  numbers  of  dissociation  and  recom¬ 
bination  reactions,  governed  by  the  local  equilibrium 
rates,  naturally  satisfies  the  microscopic  reversibility 
principle.41 

4.  Third  particle  X  is  picked  randomly  from  all  particles  in 
the  spatial  cell,  and  the  new  velocities  of  X  and  A1A2 
are  calculated  using  the  VHS/VSS  algorithm  with  the 
conservation  of  momentum  and  energy. 

For  this  algorithm  to  be  fully  defined,  we  need  to  spec¬ 
ify  the  energy  form  of  the  reaction  probability  P,  and  the  rule 
for  the  selection  of  internal  energy  levels  of  the  newly  created 
molecule.  As  discussed  below,  the  detailed  balance  is  explic¬ 
itly  used  for  the  latter,  while  the  collision  theory  for  chemical 
reactions  is  applied  to  find  Pr.  To  apply  the  collision  theory  so 


that  the  reaction  equilibrium  constant  is  captured,  it  is  neces¬ 
sary  to  present  the  dissociation  and  recombination  rates  in  an 
extended  Arrhenius  form.42 


IV.  DISSOCIATION  AND  RECOMBINATION  RATES 
OF  THE  BIAS  MODEL 

The  energy  dependent  probability  of  the  bias  dissocia¬ 
tion  model  is  simple  and  straightforward  to  implement,  but 
it  does  not  allow  its  exact  analytic  integration  to  provide 
equally  simple  temperature  dependent  reaction  rate  constant. 
This  is  mostly  due  to  the  discrete  internal  energy  modes  since 
the  integral  over  the  relative  velocities  is  complemented  by 
the  summation  over  rotational  and  vibrational  levels.  Note 
that  while  the  discrete  rotational  mode  may  be  replaced 
by  its  continuous  analog,  the  vibrational  mode  cannot  be 
simplified  this  way  due  to  large  energy  spacing  between 
low  vibrational  levels.  Also  note  that  even  fully  continu¬ 
ous  internal  modes  would  not  allow  one  to  reduce  the  bias 
model  dissociation  rate  to  the  extended  Arrhenius  form  of 
ATBe 

Even  though  the  closed  analytic  form  for  the  bias  model 
is  not  available,  the  equilibrium  reaction  rates  obtained  by  this 
model  may  be  fairly  accurately  approximated  by  a  temperature 
dependence  written  in  the  Arrhenius  form.  This  is  illustrated 
in  Fig.  1  for  the  N2+N— >N  +  N  +  N  dissociation  reaction. 
In  this  figure,  the  calculated  rate  of  the  bias  model  was  based 
on  Eq.  (3).  The  constants  Aj,  Bj,  and  E,j  of  the  Arrhenius 

—Ed 

equation,  ATBe~xr ,  which  provide  the  best  fit  to  the  calculated 
bias  model  dissociation  rates  in  the  temperature  range  between 
3000  K  and  20  000  K  are  listed  in  Table  I,  along  with  the  rec¬ 
ommended  constants  of  the  bias  model  Af/  and  Here,  the 
reaction  rate  constants  are  in  m3/( molecule/s),  and  the  ratio  of 
the  dissociation  energy  threshold  to  the  Boltzmann  constant 
Edik  is  in  K.  Figure  1  shows  that  the  Arrhenius  form  approx¬ 
imates  the  calculated  rate  very  well,  with  the  maximum  error 
in  the  considered  temperature  range  not  exceeding  5%  (and 
typically  less  than  2%).  The  approximation  is  almost  as  good 
for  the  other  dissociation  reactions  of  air  species,  also  listed  in 
Table  I. 


FIG.  1.  Reaction  rate  for  N2  +  N— >N  +  N  +  N  dissociation.  Symbols — 
calculated  rate,  line — its  Arrhenius  approximation. 
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TABLE  I.  Arrhenius  reaction  constants  matching  the  bias  model  rates  for 
different  dissociation  reactions. 


Reaction 

A 

1 

Ed/k 

N2  +  N- 

-*  N  +  N  +  N 

1.5 

2 

2.56  x  10~10 

-1 

113  400 

N2  +  N2 

-»  N  +  N  +  Nj 

1.5 

4 

8.28  x  10~u 

-1 

113  400 

02  +  O  - 

-»o  +  o  +  o 

4 

2 

1.83  x  10~10 

-1 

59370 

02  +  02 

— >  O  +  O  +  O2 

2 

4 

4.093  x  lO^11 

-1 

59370 

O?  +  Ar 

->  O  +  O  +  Ar 

0.25 

2 

1.17  x  10~u 

-1 

59370 

02  +  N2 

— >  O  +  O  +  N2 

0.4 

2 

2.50  x  10~u 

-1 

59370 

TABLE  II.  Constants  in  the  recombination  reaction  rates  for  the  bias  model. 
The  units  of  the  recombination  rate  constant  are  m6/(molecules2/s). 


Reaction 

Ar 

Br 

n2  +  n->n  +  n  +  n 

2.36  x  10~41 

-1 

N2  +  N2  -»  N  +  N  +  N2 

7.64  x  IO-42 

-1 

O2  +  0  — >0  +  0  +  0 

2.53  x  IO-43 

-0.5 

O2  +  O2  — >  O  +  O  +  O2 

5.66  x  10~44 

-0.5 

02  +  Ar-  -»  O  +  O  +  Ar 

1.62  x  10~44 

-0.5 

O2  +  N2  — *  O  +  O  +  N2 

3.46  x  IO-44 

-0.5 

The  equilibrium  constants  used  in  this  work  follow  the 
recommendations  of  Ref.  42  for  the  temperature  range  between 
3000  K  and  8000  K,  with  the  nitrogen  and  oxygen  equilibrium 
constants  K^2  and  K®2  written  as 

Kf2  =  1.084  x  1031  e“113400/rmolecule/m3,  (6) 

K °2  =  7.227  x  1032  T~0-5  e“59370/7molecule/m3.  (7) 

The  recombination  reaction  rate  constant  kr  is  then 
obtained  as  the  ratio  of  the  corresponding  dissociation  rate 
kd  to  the  equilibrium  constant  K and  may  be  presented  in  the 
form 

kr=ArTB (8) 

The  constants  Ar  and  B,  are  listed  in  Table  II.  Here,  the  dimen¬ 
sions  of  the  recombination  constant  A,.  are  m6/(molecule2/s). 
In  what  follows,  the  rates  from  Tables  I  and  II  are  used  except 
where  indicated  otherwise. 

V.  RECOMBINATION  REACTION  PROBABILITY 

As  soon  as  the  recombination  reaction  rate  is  written  in 
the  Arrhenius  form  of  Eq.  (8),  the  recombination  probability 
may  be  calculated  based  on  the  collision  theory  for  chemical 
reactions,  and  its  derivation  generally  follows  the  procedures 
outlined  in  Refs.  33  and  34.  The  main  goal  here  is  to  define  a 
recombination  probability,  Pr,  which  provides  the  number  of 


collisions  that  at  equilibrium  follows  the  recombination  rate 
listed  in  Sec.  IV,  and  thus  allows  one  to  match  the  equilibrium 
constant  at  thermochemical  equilibrium. 

The  rate  of  change  of  atomic  species  A\  due  to  recombi¬ 
nation 

Ai  +A2+X  — >  A1A2  +  X  (9) 


may  be  written  as 


dnAl 

dt 


=  nAl  «a2  nxkr. 


(10) 


where  n  is  the  number  density  of  the  corresponding  species. 
The  recombination  rate  kr  is  written  as 

K-  J  ■■■  J  c^a^fivA^fivA^tvxWv^dvAzdvx,  (11) 

(3)  (3) 

where  cr  and  oy  are  the  relative  velocity  and  reaction 
cross  section  for  the  three-body  collision,  respectively,/  is  the 
molecular  velocity  distribution  function,  and  v  is  the  veloc¬ 
ity.  The  molecular  velocities  have  X,  Y,  and  Z  components, 

and  thus  the  integral  is  nine-fold.  Note  that  the  three-body 

(3) 

recombination  reaction  cross  section  cry  has  the  dimensions 
of  m5/molecules2  (see  Ref.  43). 

Since  the  main  requirement  imposed  on  the  model  is 
matching  the  equilibrium  rate,  hereafter  we  will  assume  the 
Maxwellian  form  of  the  distribution  function  /.  It  is  also 
assumed  that  the  recombination  probability  depends  on  the  rel¬ 
ative  velocity  of  colliding  atoms,  cr,  but  not  on  the  velocity  of 
the  third  particle  X.  This  is  reasonable  since  the  internal  states 
of  the  recombining  molecules  are  selected  from  the  detailed 
balance  considerations  independent  of  the  relative  collision 
velocities  (see  Sec.  VI).  Note  also  that  at  equilibrium,  c3-3 ’  can 
be  expressed  as  a  function  of  the  ordinary  two-body  relative 

velocity  cr  and  the  masses  of  the  colliding  particles.43  In  this 

(3)  (3) 

case,  one  can  take  the  integration  over  cr  and  oy  out  of  the 
integral  over  Vx,  and  simplify  the  product  as 

OO 

[  c<;]>(Tiiy>f(vx)dvx  -  M  cr—  -  —  cr  o-tot(cr).  (12) 
J  nx  nx 

— OO 

(3) 

Here,  M  is  a  mass-dependent  constant  that  relates  cr  and 
cr,43  <Ttot(cr)  is  the  atom-atom  total  collision  cross  section, 
which  in  this  work  is  defined  by  the  VHS/VSS  intermolecular 
potential,  and  P,  is  the  recombination  probability,  which  is 
conventionally  the  ratio  of  the  reactive  to  total  collision  cross 
section,  and  in  this  case  also  includes  the  constant  M. 

Then,  substituting  Eq.  (12)  into  Eq.  (11),  and  changing 
the  variables  vAl  and  vAl  to  the  center  of  mass  velocity  ccm  and 
the  relative  velocity  cr,  as  discussed  in  Sec.  IV  C  of  Ref.  2,  for 
the  flow  at  equilibrium,  we  have 


kr 


—  Cr  tJto,{Cr)f{vAl)f{vA2)dvAidvAl 

nx 


2  (mMmAl)i 
n(kTf 


OO 

/ 


c2  e- 


(mAl+mA2)<4n 


"  dCr, 


OO 

/ 


r  r  o  _  mrcr 

- 0-wt(cr)c;  e  2 kT  dcr 

nx 


3 

2(>wAi'«a2)2 

n{kT)3 


(13) 
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where  mr,  m,\ , ,  and  ma2  are  the  reduced  mass  and  the  masses 
of  atoms  A  i  and  A2,  respectively. 

The  primary  unknown  in  this  expression  is  the  reaction 
probability,  Pr.  Let  us  now  recall  the  functional  dependence 
of  the  recombination  reaction  rate  at  equilibrium,  Eq.  (8).  Since 
the  rate  written  in  Eq.  (13)  for  equilibrium  conditions  needs  to 
be  equal  to  Ar  Tf ,  it  is  reasonable  '  '4  to  assume  the  relative 
velocity  dependent  form  of  P,  as 

Pr  =  ®cj,  (14) 


where  <F  and  VF  are  constants  to  be  found  by  equating  the  right 
hand  sides  of  Eqs.  (8)  and  (13). 

For  the  VHS/VSS  intermolecular  potentials,  the  total 
collision  cross  section  is2 


c Ttot  =  n  d 


ref' 


2  kTref\a _ 1_ 

mrcr  )  r(2  -  a) 


(15) 


Here,  E  is  the  gamma-function,  dref  and  Tref  are  the  ref¬ 
erence  diameter  and  temperature  of  the  VHS/VSS  mod¬ 
els,  respectively,  and  a  is  the  exponent  in  the  viscosity- 
temperature  dependence,  /jocT°-5+a.  For  hard  sphere  and 
Maxwell  molecules,  a  is  0  and  0.5,  respectively. 


We  also  need  to  recall  that  the  number  of  collisions  Ncou 
between  atomic  species  A  |  and  A2,  calculated  in  any  given 
collision  cell,  is  proportional  to  the  product  of  species  number 
densities  as 


Ncol,  « 


nA2 


e 

where  the  constant  e  related  to  the  number  of  unique  pairs  is 
2  when  A \  =  /L,  and  1  otherwise.  This  constant  needs  to  be 
accounted  for  when  finding  the  reaction  probability  P,  from 
Eq.  (13). 

Then,  taking  into  consideration  e  and  substituting  Eq.  (15) 
into  Eq.  (13), 


4  1 

(  Wr 

f  2kTref  \ 

l"  ndref  O 

eV/r 

\2kT)  1 

{  mr  ) 

1  r(2  -a)nx 

x  J  c]  2a+'v  e  ir  dcr.  (16) 

0 


To  integrate  Eq.  (16),  it  is  convenient  to  change  variables 

from  relative  velocity  of  colliding  atoms  cr  to  relative  collision 

2 

energy  Er,  E,  =  -Ap-,  so  that  the  integral  in  that  equation 
becomes 


/3_2q'+vF  -  mr  c r  1 

c,  e  m  dc,  - 
0 


1  l2kT 

2  \  mr 


2-a+| 

"  r 


(17) 


From  Eqs.  (16)  and  (17),  we  can  find  the  final  expression  for 
the  reaction  rate  constant, 

2  r(2~Q'+|)  2  0)  (2kT\°-5~a+^  (2kTref\a 

eyfn  T{2-a)  *  ref  nx\  mr  /  \  mr  / 

(18) 

The  constants  and  T  in  the  recombination  probability  of 
Eq.  (14)  are  found  by  equating  the  right  hand  sides  of  Eqs.  (8) 
and  (18)  and  the  exponents  in  the  temperature  dependence, 

eyfn  T(2-a)  Ar  nx  /  mr  \  '  2P 

r(2-a  +  f)r-/A4/l^i  ’  (19) 

y  =  2Br  +  2a-  1. 

The  reaction  probability  of  Eq.  (14)  is  thus  defined  and 
may  be  used  explicitly  in  the  DSMC  simulation. 

VI.  INTERNAL  ENERGY  STATES  OF  MOLECULES 
CREATED  IN  RECOMBINATION 

The  second  most  important  constraint,  after  modeling  the 
number  of  recombination  events  which  preserves  the  equi¬ 
librium  constant,  is  satisfying  the  principle  of  detailed  bal¬ 
ance  at  equilibrium.  Note  that  this  is  especially  critical  for 
the  vibrationally  favored  dissociation  model  such  as  the  bias 
model  used  here.  If  the  number  of  molecules  that  dissociate 
from  a  particular  vibrational  level  v  will  not  be  equal  to  the 


number  of  molecules  that  populate  that  level  after  recombi¬ 
nation,  then  the  vibrational  distribution  function  will  deviate 
from  its  Boltzmann  form,  and  maintaining  correct  forward  and 
reverse  reaction  rates  at  equilibrium  will  not  be  possible.  The 
main  idea  of  the  proposed  approach  is  for  each  spatial  cell  to 
record  the  rotational  and  vibrational  levels  of  molecules  that 
dissociated  in  that  cell  and  assign  them  to  molecules  created 
after  recombination. 

The  following  scheme  is  proposed  that  strictly  satisfies  the 
detailed  balance  requirement.  First,  the  rotational  and  vibra¬ 
tional  levels  of  dissociating  molecules  are  recorded  in  two 
dedicated  2D  arrays,  one  for  rotational  and  the  other  for  vibra¬ 
tional  levels  (although  those  may  be  combined  in  a  single  array 
if  needed).  The  arrays  are  cell-based,  which  means  that  each 
spatial  cell  has  its  own  list  of  rotational  and  vibrational  levels 
of  dissociated  molecules.  The  first  array  dimension  is  there¬ 
fore  dedicated  to  the  spatial  cells,  and  the  second  dimension  is 
the  memory  depth  (how  many  of  the  most  recently  dissociat¬ 
ing  molecules  will  be  remembered).  Calculations  have  shown 
that  the  memory  depth  does  not  have  to  be  large;  typically, 
storing  only  five  rotational  and  vibrational  levels  in  each  cell 
is  sufficient.  Then,  only  the  levels  of  the  last  five  dissociated 
molecules  are  stored,  the  levels  of  molecules  that  dissociated 
earlier  are  discarded.  The  arrays  are  filled  sequentially,  and 
when  the  corresponding  cell-based  running  counter  exceeds 
the  array’s  maximum  length  (the  memory  depth),  it  returns 
to  the  first  element  of  the  array,  overriding  the  values  stored 
earlier. 
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The  lists  of  stored  rotational  and  vibrational  levels  of  dis¬ 
sociated  molecules  are  used  to  populate  the  internal  states  of 
molecules  created  after  recombination.  A  cell-based  counter 
is  increased  by  one  every  time  a  pair  of  atoms  is  accepted  for 
recombination,  and  rotational  and  vibrational  levels,  J  and  v, 
which  correspond  to  the  value  of  that  counter,  are  selected  from 
the  stored  list.  Note  that  in  recombination  dominated  flows, 
there  may  be  cells  where  the  recombination  starts  before  dis¬ 
sociation,  so  that  initially  there  are  no  internal  energy  lists  of 
dissociated  molecules.  In  this  case,  one  can  take  into  consider¬ 
ation  the  functional  dependence  of  the  dissociation  probability 
of  the  bias  model,  and  consider  only  its  dominant  exponential 
term  for  the  reverse  process  of  recombination.  In  this  case,  the 
distribution  function  of  vibrational  energy  Ev  of  newly  created 
molecules  may  be  approximated  as 


The  vibrational  energy  then  sampled  as 

E„  =  ^  ln(l  +mek  -  1)),  (20) 

A, 

where  is  the  random  number  uniformly  distributed  between 
0  and  1 .  After  that,  the  vibrational  level  v  of  the  new  molecule 
is  set  to  the  maximum  value  for  which  energy  does  not  exceed 
Ev ,  and  the  rotational  level  is  set  not  to  exceed  the  remaining 
energy.  Although  the  procedure  of  Eq.  (20)  only  approximates 
the  detailed  balance,  it  is  used  only  until  the  first  dissociation 
and  thus  has  no  negative  effect  for  steady  state  flows,  and 
negligible  effect  for  transient  recombination  dominated  flows, 
as  shown  in  Sec.  VIII. 

The  after-recombination  velocities  of  particles  A i  A2  and  X 
are  computed  in  two  steps  that  provide  momentum  and  energy 
conservation.  At  the  first  step,  the  velocity  of  molecule  A 1A2  is 
set  to  the  center  of  mass  velocity  of  the  system  of  two  atoms  A  \ 
and  At,  as  dictated  by  the  conservation  of  momentum  require¬ 
ment.  The  relative  translational  energy  between  particles  A 1A9 
and  X  is  then 

eMA2-x  =  e(MA2)-x  +  ea !-a2  _  Er{J)  _  E^  _  E^  (21 ) 

Here,  E<A]Al)  x  is  the  relative  translational  energy  of  the  center 
of  mass  of  the  orbiting  pair  (A1A2)  and  the  third  particle  X, 
eA\-Ai  js  t|le  relative  translational  energy  of  the  two  atoms, 
and  Er(J )  and  Ev(v)  are  the  rotational  and  vibrational  energies 
of  levels  J  and  v,  respectively.  Note  that  for  some  J,  v,  and 
e(AiA2)-x^  tjje  vaiue  0f  eA\A2-x  may  |1C  |ess  tjlan  0  ]n  this 

case,  another  particle  X  needs  to  be  randomly  selected  from 
all  particles  in  the  cell,  while  the  selected  levels  J  and  V  should 
not  be  changed.  Alternatively,  all  three  particles  At,  A2,  and  X 
may  be  re-selected. 

At  the  second  step,  the  final  velocities  of  A 1  As  and  X  are 
calculated  using  the  conventional  VHS  or  VSS  algorithm2  for 
the  relative  translational  energy  Ey1  2~x . 

VII.  VERIFICATION  OF  THE  RECOMBINATION  MODEL 

A.  Detailed  balance  at  thermal  and  chemical 
equilibrium 

For  any  dissociation  model,  the  verification  of  the  cor¬ 
responding  recombination  model  is  fairly  straightforward: 


one  needs  (i)  to  verify  that  the  model  is  capable  of  main¬ 
taining  thermochemical  equilibrium  with  the  reactant  mole 
fractions  that  correspond  to  the  given  equilibrium  constant, 
in  a  wide  range  of  temperatures  and  pressures,  and  (ii) 
check  that  the  system  that  was  initially  at  non-equilibrium 
converges  to  the  correct  equilibrium  condition.  Equilibrium 
here  needs  to  be  both  microscopic  (Maxwellian  veloc¬ 
ity  distributions  and  Boltzmann  internal  energy  popula¬ 
tions)  and  macroscopic  (equal  mode  temperatures  and  the 
ratio  of  species  mole  fractions  matching  the  equilibrium 
constant). 

With  this  in  mind,  the  first  set  of  tests  was  conducted 
for  a  system  initially  at  full  thermal  and  chemical  equi¬ 
librium,  with  the  concentrations  of  molecular  and  atomic 
species  set  to  their  equilibrium  value  at  a  given  tempera¬ 
ture  and  with  the  Maxwell-Boltzmann  distribution  of  parti¬ 
cle  velocities  and  internal  energies  at  that  temperature.  The 
SMILE  DSMC  code44  was  used  in  all  computations,  extended 
to  include  the  proposed  bias  model.  The  calculations  were 
conducted  for  N2-N,  O2-O,  and  CL-O-Ar  systems,  with 
dissociation-recombination  reactions  N2+M<->N  +  N  +  M  and 
02+M«->0  +  0  +  M,  where  M  is  the  third  particle  that 
can  be  of  any  species  of  the  gas  mixture.  Modeling  results 
in  an  adiabatic  reservoir  with  the  specular  gas-surface  inter¬ 
action,  gas  temperature  varying  from  3000  K  to  14  000  K, 
and  pressures  from  0.01  to  50  atm  have  indicated  that  the 
bias  dissociation  and  recombination  model  does  not  disturb 
an  equilibrium  system.  The  observed  deviation  from  equilib¬ 
rium  constant  was  less  than  0.3%.  That  deviation,  aside  from 
the  statistical  scatter  on  the  order  of  0.1%,  was  also  due  to  the 
1  %-2%  difference  between  the  actual  dissociation  rates  for  the 
bias  model  and  their  approximation  written  in  the  Arrhenius 
form. 

An  example  of  thermal  relaxation  from  the  equilibrium 
state  is  shown  in  Fig.  2(a),  where  the  translational  and  internal 
temperatures  (Ttm,  T rol ,  and  7’,,.,/,)  and  the  species  mole  frac¬ 
tions  (X[0]  and  XfCL])  are  shown  as  functions  of  time.  The 
system  is  O2-O,  the  initial  temperature  is  5000  K,  and  the  num¬ 
ber  density  is  2  x  1025  molecule/m3,  with  the  mole  fractions  of 
81.4%  0-18.6%  O2.  The  total  number  of  simulated  molecules 
in  the  system  was  approximately  ten  million,  the  total  relax¬ 
ation  time  amounted  to  about  10  000  mean  collision  times, 
and  the  number  of  dissociation  and  recombination  reactions 
during  the  relaxation  time  exceeded  two  million.  The  figure 
clearly  shows  that  the  system  does  not  deviate  from  its  initial 
state.  The  distribution  functions  are  equilibrium,  as  illustrated 
in  Fig.  2(b)  for  the  vibrational  (VDF)  and  rotational  (RDF) 
distribution  functions. 

B.  Relaxation  from  a  chemically  non-equilibrium  state 

The  second  set  of  tests,  also  conducted  for  the  above 
three  gas  systems,  was  aimed  at  the  verification  of  the  abil¬ 
ity  of  an  initially  non-equilibrium  system  to  reach  equilibrium. 
With  varying  initial  gas  pressures,  temperatures,  and  composi¬ 
tions,  the  system  was  tracked  until  the  steady  state  is  reached, 
and  then  the  velocity  and  energy  distribution  functions  and 
species  mole  fraction  were  compared  to  the  corresponding 
equilibrium  values.  An  example  of  the  temporal  relaxation 
of  a  non-equilibrium  system  is  presented  in  Fig.  3(a)  for  a 
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Rotational  Level 


(a)  (b) 


FIG.  2.  Equilibrium  relaxation  in  O2- 
O  bath:  macroparameters  (left)  and 
internal  energy  distribution  functions 
(right). 


Rotational  Level 


(a)  (b) 


FIG.  3.  Collisional  relaxation  toward 
equilibrium  state  (left)  and  internal 
energy  distribution  functions  at  6  ns 
(right).  N2-N  thermal  bath. 


N2-N  gas  mixture.  Initially,  the  gas  is  100%  N2  with  a  tem¬ 
perature  of  14  000  K  in  all  modes  and  a  number  density  of 
1027  molecule/m3.  The  collisional  relaxation  results  in  a  pref¬ 
erential  dissociation  in  the  first  nanosecond,  followed  by  slow 
relaxation  toward  thermochemical  equilibrium.  The  equilib¬ 
rium  is  reached  at  about  6  ns.  The  internal  energy  level  popu¬ 
lations  at  that  time  become  fully  equilibrium,  as  illustrated  in 
Fig.  3(b). 


VIII.  RECOMBINATION  MODEL  VALIDATION 

Chemistry  model  validation  is  traditionally  difficult  for 
kinetic  methods,  mostly  due  to  limited  availability  of  exper¬ 
imental  data  in  near-continuum  flow  regimes  where  kinetic, 
non-equilibrium  effects  would  be  pronounced.  Validation  is 
even  more  challenging  for  recombination  models  since  most 
shock  tube  data,  being  the  primary  scope  of  DSMC  stud¬ 
ies  dealing  with  dissociation  and  exchange  reactions,  are 
dissociation-dominated,  and  the  impact  of  recombination  is 
relatively  minor  and  generally  too  small  for  detailed  validation. 
In  this  work,  we  have  chosen  ozone  pyrolysis47  and  photol¬ 
ysis46  experiments  where  the  ozone  molecules  are  quickly 
transformed  to  oxygen  atoms,  which  then  recombine  through 
collisions  with  the  surrounding  gas  to  produce  molecular  oxy¬ 
gen.  In  the  first  case,47  pyrolysis  of  ozone  was  used  to  produce 
an  excess  of  O  atoms  in  the  shock  tube,  and  oxygen  col¬ 
lisional  recombination  then  proceeded  in  a  thermal  bath  of 


argon  heated  to  temperatures  between  2000  K  and  3000  K.  In 
the  second  case,46  experiments  were  conducted  at  room  tem¬ 
peratures  in  500  Torr  nitrogen,  with  metered  portions  of  the 
flow  picking  up  0.25-0.75  Torr  ozone,  which  was  then  almost 
fully  dissociated  by  the  excimer  laser  beam,  and  the  resulting 
oxygen  atoms  recombined  through  three-body  collisions  with 
nitrogen. 

Numerical  simulation  of  gas  conditions  of  the  shock  tube 
experiments4  is  conducted  here  for  a  0.5%  0-99.5%  Ar  mix¬ 
ture  initially  at  40  cm  Hg  and  2400  K.  Note  that  the  argon 
temperature  increased  somewhat  due  to  heat  release  in  recom¬ 
bination,  but  that  increase  was  only  about  20  K.  As  the  ref¬ 
erence  point,  the  first  computation  was  conducted  with  the 
dissociation  reactions  turned  off.  In  this  case,  only  the  recom¬ 
bination  process,  driven  by  the  recombination  rate  at  that  tem¬ 
perature,  impacts  the  oxygen  mole  fractions.  This  makes  the 
result  comparable  to  the  recombination-only  DSMC  calcula¬ 
tion  of  Ref.  8.  For  full  compatibility,  these  recombination-only 
computations  have  used  the  recommended3  oxygen-argon 
recombination  rate  with  Ar  =  1.1  X  10~44  m6/(molecule2/s). 
The  comparison  with  Ref.  8  is  given  in  Fig.  4(a),  where  the 
present  results  are  labeled  “Modeling,  R,”  and  those  of  Ref. 
8,  “Koura  (1994).”  The  plot  shows  the  temporal  relaxation 
of  the  ratio  of  the  initial  number  density  of  oxygen  atoms, 
no( 0),  to  the  instantaneous  time-dependent  oxygen  atom  den¬ 
sity  «o(f).  Clearly,  there  is  good  agreement,  and  a  small 
difference  is  attributed  to  the  statistical  scatter  inherent  in 
DSMC. 
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FIG.  4.  Temporal  relaxation  of  O  pop¬ 
ulation  in  Ar  (left)  and  No  (right)  baths. 


Comparison  with  the  measurements  of  Ref.  45  for  a 
recombining-dissociating  gas  is  also  shown  in  Fig.  4(a).  In 
this  case,  the  reaction  rates  from  Tables  I  and  II  were  used. 
The  experimental  data  are  labeled  “Wray  (1963).”  Note  that 
the  experimental  points  are  much  lower  than  the  results  for 
the  recombination-only  flow.  The  much  slower  recombination 
process  in  the  experiment  is  attributed  primarily  to  the  disso¬ 
ciation  of  newly  created  atoms.  Even  though  the  translation- 
rotation  temperature  of  the  gas  is  relatively  low,  2400  K,  it  is 
still  sufficient  to  dissociate  oxygen  molecules  from  the  upper 
vibrational  levels.  The  population  of  such  molecules  is  dispro- 
portionally  high  since  the  recombination  mostly  leads  to  the 
formation  of  vibrationally  excited  molecules. 

The  present  baseline  simulation  is  denoted  “Modeling, 
Z„(Tr).”  In  the  baseline  setup,  the  vibration-translation  (VT) 
energy  transfer  in  Oo-Ar  collisions  was  modeled  using  the 
conventional  LB  approach  with  the  translational  temperature 
dependent  vibrational  relaxation  number  ZV=ZV(TV )  deter¬ 
mined  by  the  Millikan- White  empirical  correlation47  with 
Park’s  high  temperature  correction.44  As  clearly  seen  in  Fig. 
4(a),  the  mole  fraction  ratio  for  the  baseline  bias  model  is  much 
lower  than  that  in  the  experiment.  The  difference  is  mostly  due 
to  the  inapplicability  of  the  standard  LB  procedure  in  this  case. 
The  VT  energy  transfer  rate  computed  using  the  translational 
temperature  of  298  K  is  negligible  compared  to  the  recom¬ 
bination  and  dissociation  rates,  and  there  is  no  de-excitation 
of  oxygen  molecules  that  populate  high  vibrational  levels 
after  recombination.  As  the  result,  the  recombined  molecules 
quickly  dissociate,  and  the  population  of  CB  becomes  almost 
frozen  after  40  /-is. 

In  order  to  overcome  this  obvious  issue,  a  simple  correc¬ 
tion  to  the  LB  approach  has  been  used,  where  the  Millikan- 
White  expression  is  now  based  on  the  vibrational  and  not 
translational  temperature.  The  result,  shown  as  “Modeling, 
Z„(T„),”  indicates  significant  improvement  as  compared  to 
the  standard  LB  approach.  Still,  the  computed  rate  of  molec¬ 
ular  oxygen  formation  visibly  underpredicts  the  measured 
rate.  In  an  attempt  to  alleviate  the  problem,  a  state-to-state 
3D  Forced  Harmonic  Oscillator  (FHO)  model  9  is  also  used 
for  the  VT  energy  transfer.  The  results  obtained  using  this 
model  are  denoted  “Modeling,  FHO.”  The  3D  FHO  model 
has  been  shown  to  match  trajectory  calculations  for  all 
vibrational  levels49  and  is  believed  to  provide  adequate 


accuracy  to  the  VT  energy  exchange  of  diatomic  molecules. 
As  shown  in  Fig.  4(a),  the  use  of  that  model  results  in  good 
agreement  with  the  experimental  data.  The  computed  O2 
formation  rate  is  approximately  20%  lower  than  the  mea¬ 
sured  rate.  There  may  be  several  reasons  for  the  difference. 
First,  the  bias  dissociation  model  overpredicts  its  Arrhe¬ 
nius  fit  for  Oi-Ar  by  about  10%  at  temperatures  below 
3000  K.  Second,  the  equilibrium  constant42  becomes  less 
accurate — and  generally  too  high — at  temperatures  below 
3000  K.  Finally,  the  recombination  model  proposed  in  this 
work,  while  accurate  in  providing  detailed  balance  near  equi¬ 
librium,  may  be  less  accurate  for  highly  nonequilibrium  flows, 
where  it  is  not  clear  how  well  the  assignment  of  inter¬ 
nal  energies  of  dissociating  molecules  to  molecules  created 
after  the  recombination  reproduces  the  actual  reaction  cross 
sections. 

The  second  series  of  validation  computations  have  been 
conducted  for  the  conditions  that  reproduce  those  of  the  Ref.  46 
experiments.  Atomic  oxygen,  initially  at  a  partial  pressure  of 
0.57  Torr,  is  modeled  in  504  Torr  N2  at  a  room  temperature  of 
298  K  (this  temperature  increased  with  time  to  302  K  because 
of  the  recombination).  The  major  difficulty  of  comparing  to 
the  data46  is  the  largely  unknown  oxygen  dissociation  rate 
at  such  a  low  temperature.  In  this  case,  it  is  reasonable  to 
use  the  recommended46  oxygen  in  the  nitrogen  recombina¬ 
tion  rate  of  kr  -  1 .9  x  10-45  m6/(molecule2/s)  and  calculate  the 
corresponding  dissociation  rate  using  the  300  K  equilibrium 
constant  of  5  x  10-56  molecule/m3.  In  this  case,  the  parameter 
Ad  of  the  bias  dissociation  model  is  approximately  0.05,  and 
not  0.4  valid  for  higher  temperatures. 

Comparison  of  the  computed  atomic  oxygen  mole  frac¬ 
tion  with  the  experiment46  is  therefore  given  in  Fig.  4(b)  for 
three  different  cases:  the  recombination-only  case  (“Modeling, 
R”),  the  recombination-dissociation  case  with  Ad  =0.4  (“Mod¬ 
eling,  RDV,  A  =  0.4”),  and  the  recombination-dissociation 
case  with  Ac/  =  0.05  (“Modeling,  RDV,  A  =  0.05”).  The  3D 
FHO  model50  was  used  for  the  VT  energy  transfer.  As 
expected,  the  recombination-only  relaxation  follows  closely 
the  analytical  slope  of  the  fitted  experimental  rate  of  kr 
=  1.9  x  10-45  m6/(molecule2/s).  The  dissociation  predictably 
decreases  that  slope.  Note  that  the  experimental  fit  proposed 
in  Ref.  46  does  not  take  into  account  the  possibility  of  the 
concurrent  dissociation  from  high  vibrational  levels,  and 
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thus  the  difference  between  the  experimental  fit  and  the  full 
dissociation-recombination  model  could  be  expected.  That  dif¬ 
ference  is  small,  however,  for  the  realistic  dissociation  rate  of 
Aj  =  0.05.  The  use  of  the  high-temperature  dissociation  rate 
results  in  the  underprediction  of  the  experimental  slope  by  over 
30%. 

IX.  CONCLUSIONS 

Kinetic  simulations  of  high  temperature  flows  often  need 
fast  and  reliable  models  that  capture  the  key  physics  and  that 
can  be  used  on  reactions  for  which  detailed  cross-sectional 
information  is  not  available.  The  bias  dissociation  model  is 
unique  in  needing  no  a  priori  assumption  of  an  Arrhenius 
rate  coefficient  and  essentially  no  adjustable  parameters  (its 
pre-exponential  factor  is  on  the  order  of  unity).  It  reproduces 
vibration-dissociation  coupling,  equilibrium,  and  nonequilib¬ 
rium  behavior  as  a  natural  consequence  of  its  physics-based 
cross-sectional  function  and  just  a  simple  chemistry-based 
assessment  of  whether  the  reaction  will  be  strongly  favored 
or  not.  One  disadvantage  of  the  model  used  to  be  the  lack  of  a 
compatible  recombination  model. 

To  fill  this  gap,  a  recombination  model  is  developed  in  this 
work  suitable  for  the  direct  simulation  Monte  Carlo  method 
and  compatible  with  the  bias  dissociation  model.  The  model 
captures  the  temperature  dependence  of  the  equilibrium  reac¬ 
tion  constant  and  satisfies  the  detailed  balance  requirement  at 
equilibrium.  The  recombination  probability  is  derived  using 
the  collision  theory  for  chemical  reactions,  and  the  rotational 
and  vibrational  modes  of  newly  created  molecules  are  pop¬ 
ulated  with  the  corresponding  levels  of  most  recently  disso¬ 
ciated  molecules.  The  model  is  straightforward  to  implement 
and  was  shown  to  provide  reasonable  agreement  with  results 
of  previous  numerical  simulations  and  experiments  on  atomic 
oxygen  recombination. 

The  approach  to  modeling  the  recombination  process  is 
general  enough  to  be  used  for  many  other  dissociation  models. 
While  the  present  implementation  assumes  the  independence 
of  the  internal  energy  states  of  molecule  created  after  recombi¬ 
nation  on  the  chemical  species  of  the  third  particle,  that  depen¬ 
dence  may  easily  be  incorporated.  It  would  only  require  the 
extension  of  the  array  that  stores  rotational  and  vibrational  lev¬ 
els  of  dissociated  molecules  to  include  the  information  about 
the  third-particle  species.  Note  also  that  the  present  imple¬ 
mentation  considers  only  the  recombination  of  two  atoms  into 
a  diatomic  molecule.  While  the  extension  to  the  dissociation 
and  recombination  of  polyatomic  molecules  is  possible,  it  may 
be  difficult  due  to  significantly  higher  memory  requirements 
for  the  array  of  internal  states,  as  well  as  the  need  to  take 
into  consideration  the  internal  states  of  the  molecular  products 
of  dissociating  polyatomic  molecules  and  molecular  reactants 
that  participate  in  recombination. 
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